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Abstract Bayesian tip dating is a recently developed method to estimate divergence times and 
evolutionary rates. It overcomes several drawbacks in traditional stepwise approach. However, it 
also requires more knowledge about statistics. This paper hierarchically explains the theory and 
computation in the Bayesian tip-dating approach, and divides the whole process into prior for the 
divergence times, prior for the evolutionary rates, model for the character changes and Markov 
chain Monte Carlo algorithm, which are key components in this method. The aim is to provide a 
general guidance for paleontologists in empirical data analyses. 

Key words Bayesian tip dating, fossilized birth-death process, relaxed clock, Mk model, 
MrBayes 


Citation Zhang C, in press. Using Bayesian tip-dating method to estimate divergence times and 
evolutionary rates. Vertebrata PalAsiatica. DOI: 10.19615/j.cnki.2096-9899.210516 


贝 叶 斯 支 端 定年 法 推断 分 异 时 间 和 演化 速率 


张驰” 
(1 中 国 科 学 院 古 疹 椎 动物 与 古人 类 研究 所 ， 中 国 科 学 院 疹 椎 动物 演化 与 人 类 起 源 重点 实验 室 北京 100044) 
2 中 国 科学 院 生 物 演化 与 环境 卓越 创新 中 心 北京 100044) 


摘要 : 贝 叶 斯 文 端 定 年 法 是 近 些 年 开发 的 推断 类 和 群 分 异 时 间 和 演化 速率 的 方法 。 它 克服 了 
传统 分 步 计 算 的 缺陷 ， 但 涉及 到 的 统计 学 知识 也 更 多 。 本 文 从 贝 叶 斯 统计 计算 的 角度 分 层 
剖析 了 支 端 定年 法 的 原理 和 计算 过 程 ， 按 照 分 异 时 间 的 先 验 分 布 、 演 化 速率 的 先 验 分 布 、 
寺 征 状 态 变 化 的 模型 和 马 氏 链 蒙 特 卡 罗 算 法 几 个 部 分 ， 叙述 并 讨论 了 定年 计算 中 的 主要 模 
型 和 算法 。 旨 在 一 定 程度 上 为 古生物 学 家 分 析 实 际 数据 提供 参考 。 
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推 新 类 和 群 的 系统 关系 和 分 异 时 间 是 分 文系 统 学 分 析 的 基础 。 如 何 合理 利用 化 石 的 形 
态 和 地 质 年 代 的 数据 来 完成 此 类 推断 一 直 是 一 个 环 手 的 问题 。 传 统 的 推 朵 过 程 往往 采用 
分 步 的 策略 。 首 先 单 独 从 形态 数据 出 发 ， 利 用 简约 法 构建 类 群 的 系统 关系 。 这 个 关系 只 
有 拓扑 结构 而 没有 时 间 的 信息 ， 而 且 只 代表 给 定时 间 和 搜索 条 件 下 最 优 的 结果 ( 即 最 简约 
的 树 或 其 合意 树 )。 然 后 固定 这 个 拓扑 结构 ， 从 化 石 年 代 出 发 ， 利 用 最 小 校长 法 (Laurin， 
2004) 或 等 距 枢 长 法 (Brusatte et al., 2008) 确 定 树 中 内 部 节点 的 分 异 时 间 。 最 小 枝 长 法 是 把 当 
前 节点 的 时 间 往 前 推 一 个 百 万 年 作为 祖先 节点 的 时 间 ; 而 等 踢 枝 长 法 则 是 把 祖先 节点 到 
后 代 节 点 的 时 间 的 中 点 最 为 当前 节点 的 时 间 。 至 于 演化 速率 ， 可 以 通过 祖先 状态 重建 和 
每 个 树 校 上 特征 变化 的 次 数 ， 结 合 上 一 步 推断 的 时 间 来 估计 。 这 种 分 步 计 算 的 方式 比较 
直观 ， 因 此 被 用 于 实际 数据 分 析 (Wang and Lloyd, 2016)。 然 而 ， 该 策略 存在 诸多 缺陷 。 首 
先 ， 它 每 一 步 都 忽略 了 推断 中 的 不 确定 性 ， 包 括 树 的 拓扑 结构 、 分 异 时 间 和 祖先 特征 状 
aS; 其 次 ， 每 一 步 都 只 利用 了 一 部 分 数据 信息 ， 如 建树 时 只 用 到 形态 特征 ， 定 年 时 只 用 
到 化 石 年 代 ; 再 次 ， 定 年 的 方式 很 主观 ， 对 化 石 数量 的 增 减 很 敏感 ， 且 不 适用 于 现 生 类 
群 ; 最 后 ， 整 个 过 程 缺乏 一 个 严谨 的 统计 学 框架 ， 无 法 对 不 同 的 模型 假设 进行 检验 。 

近 些 年 开发 的 贝 叶 斯 支 端 定年 法 (Bayesian tip dating) (Ronquist et al., 2012; 
Gavryushkina et al., 2014; Zhang et al., 2016) 很 好 地 克服 了 上 述 问 题 。 贝 叶 斯 支 端 定 年 法 
把 化 石 形态 和 年 代数 据 整 合 在 一 次 完整 的 计算 过 程 中 ， 能 够 尽 可 能 地 利用 数据 信息 ， 同 
时 考虑 了 树 的 拓扑 结构 、 分 异 时 间 、 演 化 速率 以 及 化 石 年 代 的 不 确定 性 。 该 方法 通过 统 
计 模 型 来 描述 特征 的 演化 、 类 群 的 生 灭 以 及 化 石 的 采样 等 过 程 ， 并 借助 相对 成 熟 的 贝 叶 
斯 统计 框架 和 计算 方法 来 进行 参数 估计 和 模型 选择 。 但 该 方法 相对 来 说 比较 复杂 ,需要 
较 多 的 统计 学 知识 ， 对 古生物 学 家 来 说 往往 难于 理解 和 上 手 ， 而 且 人 鲜 有 系统 阐述 支 端 定 
年 法 计算 过 程 和 参数 意义 的 文献 (Gavryushkina and Zhang, 2020)。 本 文 逐 层 训 析 文 端 定 
年 法 的 计算 过 程 ， 解 释 其 中 用 到 的 重要 模型 和 参数 意义 ， 则 在 一 定 程度 上 为 古生物 学 家 
分 析 实 际 数据 提供 参考 。 

本 文 首先 介绍 描述 时 间 树 的 石化 生 灭 过 程 (fossilized birth-death process) 模 型 (Stadler, 
2010), 然后 介绍 描述 特征 演化 速率 的 宽松 形态 钟 (relaxed clock) 模 型 ， 接 着 介绍 描述 特征 
状态 变化 的 Mk 模型 (Lewis, 2001), 再 通过 贝 叶 斯 公式 把 上 述 模 型 联系 起 来 ， 最 后 介绍 估 
计 参 数 后 验 分 布 的 马 氏 链 蒙特 卡 罗 (Markov chain Monte Carlo, MCMC) 算 法 。 附 录 提 供 
了 计算 中 生 代 乌 类 数据 (Zhang and Wang, 2019) 的 MrBayes 命 令 。 


2 分 异 时 间 


时 间 树 (timetree) 代 表 类 群 的 系统 关系 和 分 异 时 间 ， 它 的 概率 分 布 可 以 通过 石化 生 
灭 过 程 (Stadler 2010) 来 给 出 。 该 过 程 描 述 了 从 这 些 类 群 的 最 近 共 同 祖先 ( 树 根 ) 开 始 ， 
分 异 、 灭 绝 、 采 集 化 石和 采样 现 生 类 群 这 一 系列 事件 的 发 生 ， 并 对 应 于 一 棵 完整 树 ( 图 
1A)。 但 是 实际 数据 分 析 中 无 法 推 呆 这 个 完整 树 ， 只 能 推断 和 样本 相关 的 部 分 ， 即 样本 
树 (图 1B)。 记 每 个 树枝 的 分 异 速率 (或 叫 成 种 率 ) 为 4 灭绝 速率 为 沿 每 个 树枝 的 化 石 采 


202105.00071v1 


chinaXiv 


ChinaXiv 合 作 期 刊 


Zhang - Using Bayesian tip-dating method 3 


样 速率 为 y, 现 生 类 群 的 采样 概率 (或 采样 比例 ) 为 p, 通过 建立 并 求解 一 系列 常 微分 方程 ， 
可 以 得 到 给 定 4, u, y, p 时 ， 样 本 时 间 树 T= {c, 丰 的 概率 分 布 ， 记 为 P(T | 4, u, yw, p), 其 中 zt 


代表 拓扑 结构 ，t 代 表 以 百 万 年 为 单位 的 分 异 时 间 。 


A 


B 


图 1 石化 生 灭 过 程 可 能 得 到 的 时 间 树 
A. 生 灭 过 程 对 应 的 完整 树 ; 
B. 只 保留 和 样本 相关 的 分 支 得 到 的 样本 树 ， 包 含 两 个 现 生 类 群 ( 红 点 ) 和 两 个 化 石 ( 蓝 点 ) 
Fig. 1 Example timetree generated from the fossilized birth-death process 
A. compete tree by keeping all branches; B. sampled tree by only keeping branches leading to two 
extant taxa (red dots) and two fossils (blue dots) 


在 MrBayes 软 件 中 ， 该 生 灭 过 程 以 树 根 时 间 
为 起 始 条 件 ( 图 2), 计算 时 要 指定 t 的 先 验 分 布 。 
通常 这 个 先 验 比较 宽泛 (最 大 范围 从 0 到 无 穷 ), 不 
过 一 般 能 从 人 研究 的 类 群 预 估 一 个 更 精确 的 范围 ， 
例如 ， 其 下 界 不 会 早 于 最 古老 化 石 的 年 代 。 化 石 
的 时 间 可 以 固定 为 具体 的 数值 ( 百 万 年 前 ), 也 可 
以 用 一 个 均匀 分 布 给 出 时 间 的 上 下 界 。 在 计算 时 
还 需要 提供 现 生 类 群 大 致 的 采样 比例 (p)。 现 生 
类 群 可 以 有 两 种 采样 策略 ， 一 种 是 均匀 随机 采样 
(random), 另 一 种 是 多 样 化 采样 (diversity) (Zhang 
et al., 2016), 可 以 根据 实际 数据 的 情况 自行 选择 ， 
后 者 可 能 更 符合 高 阶 元 类 群 的 采样 模式 (如 每 个 
科 只 取 一 个 代表 的 属 或 者 每 个 属 取 一 个 代表 物 
种 )。 对 于 分 异 、 灭 绝 和 化 石 采 样 速率 ， 程 序 为 
了 设置 先 验 方便 ， 重 新 参数 化 为 4 = 4 一 1,v =/ 
A,s=yl (4tyw)。q 的 默认 先 验 为 指数 分 布 (范围 从 
0 到 无 穷 ), v 和 s 的 默认 先 验 为 均匀 分 布 (范围 从 0 到 
1, 更 一 般 地 为 贝塔 分 布 )。 这 样 ， 时 间 树 包括 分 
异 时 间 等 参数 的 先 验 分 布 基本 就 确定 了 。 

需 特 别提 到 的 是 ， 有 些 数据 只 包含 了 化 石 
而 没有 现 生 类 群 。 这 时 一 般 假 设 生 灭 过 程 在 未 
到 达 现 今 时 间 点 所 有 类 和 群 就 都 灭绝 了 。 因 此 现 


Si :1 So :1 

图 2 用 于 各 个 概率 分 布 公式 中 参数 的 示例 
树 上 化 石 Fi 和 FEF: 的 特征 状态 为 0, 现 生 类 群 
Si 和 8S: 的 特征 状态 为 1, 内 部 节点 的 特征 状 

态 用 xo, xu zx 表示 。 
化 石 的 年 代 分 别 为 = 100 Ma, t; = 50 Ma, 
树 根 的 时 间 为 t, 其 他 分 异 时 间 为 6 和 4。 
记 t= {h, ty b, tp 1s}。 各 个 校 上 的 特征 演化 

速率 为 r= {rp Fa 73, Fa Fs} 


Fig. 2 Example parameters and symbols 


used in the probability distributions 
The character states for fossils F and F, are 
0, for extant taxa S, and S, are 1, for internal 
nodes are xo, Xi, X2. The ages of fossils F, and 
F, are t; = 100 Ma and t; = 50 Ma 
The root age is ¢, and the remaining 
divergence times are t, and t,. Denote t = {t,, 
ty, b, ta, ts}. The evolutionary rates on the 
branches are r = {71, 7, 1°35 Fa, Ts} 
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生 类 群 不 论 采 用 何 种 采样 策略 和 比例 (程序 默认 为 1.0), 都 没有 样本 被 采集 到 。 不 过 由 于 
MrBayes 软 件 的 限制 ， 最 年 轻 的 样本 总 被 显示 为 现 生 类 群 (时 间 为 0), 因此 需要 把 整 棵 树 
的 时 间 轴 进行 相应 的 平移 ， 或 者 说 分 异 时 间 要 加 上 年 代 最 近 的 化 石 的 年 代 。 这 只 是 结果 
显示 的 问题 而 不 是 程序 计算 的 错误 (后 续 版 本 将 修复 这 个 问题 )。 

除了 石化 生 灭 过 程 之 外 ，MrBayes 还 为 时 间 树 提供 了 一 个 均匀 分 布 的 先 验 (Ronquist 
et al., 2012)。 它 没有 生火 和 采样 这 些 参 数 ， 只 依赖 于 树 根 时 间 4, 因此 只 需 设 置 的 先 验 
分 布 。 均 匀 分 布 时 常 被 认为 是 没有 信息 的 先 验 ,但 它 实际 上 往往 带 有 很 强 的 信息 ， 在 定 
年 这 个 问题 上 ， 会 影响 对 分 异 时 间 的 估计 。 而 石化 生 灭 过 程 看 似 参数 很 多 ， 但 实际 上 其 
设置 可 以 很 灵活 。 例 如 ， 分 异 、 灭 绝 和 化 石 采 样 速 率 都 可 以 随时 间 变 化 ， 在 不 同 的 时 间 
段 内 各 自 独立 (Gavryushkina et al., 2014; Zhang et al., 2016)。 这 可 能 更 符合 实际 生物 学 过 
程 ， 同 时 还 能 推测 净 成 种 速率 和 化 石 采样 速率 随时 间 的 变化 。 


3 ”演化 速率 


形态 特征 的 演化 速率 是 指 每 百 万 年 每 个 特征 期 望 的 变化 次 数 。 对 于 给 定 的 一 段 时 
间 ， 演 化 速率 越 快 ， 则 特征 最 终 期 望 的 变化 次 数 越 多 。 一 般 给 每 个 树枝 一 个 演化 速率 参 
数 ， 记 为 r (图 2)。 时 钟 模 型 应 用 于 形态 数据 时 被 称 为 形态 钟 模型 ， 类 比 于 用 在 分 子 数据 
时 的 分 子 钟 模型 。 严 格 钟 (strict clock) 模 型 假设 演化 速率 在 各 个 树 桔 上 都 相同 ， 这 通常 不 
适用 于 形态 数据 ,实际 分 析 时 往往 需 使 用 宽松 钟 (relaxed clock) 模 型 。 宽 松 钟 模型 可 以 分 
为 两 类 ， 一 类 为 独立 速率 ， 男 一 类 为 自 相关 (autocorrelated) 速 率 ， 区 别 在 于 r 的 概率 分 布 
P(r) 不 同 。 

独立 速率 模型 假设 校 上 的 演化 速率 彼此 独立 ,它们 都 服从 均值 相同 的 某 个 概率 分 
布 。 常 用 的 概率 分 布 包括 伽 马 分 布 (Lepage et al., 2007) 和 对 数 正 态 分 布 (Drummond et al., 
2006)。 分布 的 均值 也 被 称 为 基准 速率 (base rate), 反映 平均 的 演化 速率 大 小 。 分 布 的 方差 
则 反映 演化 速率 在 树枝 之 间 变 化 的 剧烈 程度 : 方差 较 小 时 各 个 速率 相差 不 大 ， 这 意味 着 
演化 速率 在 整 棵 树 上 没有 明显 的 差异 ; 而 方差 越 大 ， 不 同 树 校 上 速率 的 差异 越 明 显 。 

自 相 关 速 率 模型 假设 后 代 树 校 上 的 演化 速率 依赖 于 临近 祖先 那 校 上 的 速率 (例如 产 
和 zs 都 依赖 ,7 和 都 依赖 r,)。 当 前 树 校 上 的 速率 一 般 假设 服从 对 数 正 态 分 布 (Kishino et 
al., 2001; Thorne and Kishino, 2002), 其 均值 为 临近 祖先 节点 的 速率 。 同 理 ， 分 布 的 方差 
也 反映 演化 速率 在 树枝 之 间 变 化 的 剧烈 程度 。 

这 两 类 速率 模型 往往 对 分 异 时间 的 估计 也 有 影响 ， 这 主要 是 因为 自 相 关 速 率 模型 
会 倾向 于 速率 的 变化 是 渐进 的 ， 而 独立 速率 模型 没有 这 种 限制 ， 会 更 适应 临近 树枝 间 速 
率 变化 比较 剧烈 的 情况 。 对 化 石 形态 数据 来 说 ， 独 立 速率 模型 可 能 更 适用 。 
默认 的 情况 下 ， 形态 数据 矩阵 中 所 有 特征 都 共享 每 个 枝 上 的 演化 速率 ， 因 此 ， 这 
个 速率 代表 的 是 所 有 特征 的 平均 情况 。 如 果 需 要 考虑 不 同 特征 演化 速率 的 异 质 性 ， 就 
需要 对 特征 进行 分 区 。 一 般 可 以 按照 不 同 特征 类 型 或 不 同 身体 部 位 或 功能 来 分 ， 每 个 
分 区 内 的 特征 共享 一 组 演化 速率 ， 而 分 区 之 间 特 征 演化 速率 的 模式 是 独立 的 ， 这 样 就 
可 以 推断 不 同 部 位 或 功能 相关 特征 随时 间 会 发 生 怎样 的 变化 (Lee, 2016; Zhang and Wang, 
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2019)。 需 注意 的 是 ， 分 区 越 多 ， 每 个 分 区 内 的 特征 数量 就 越 少 ， 因 此 能 够 佑 计 演 化 速 
率 参数 的 信息 就 越 少 ， 会 造成 方差 很 大 甚至 参数 个 数 超过 特征 数量 导致 无 法 进行 参数 估 
计 。 因 此 ， 要 在 考虑 演化 速率 异 质 性 和 分 区 数量 之 间 做 一 个 权衡 。 


4 ”特征 状态 变化 


有 了 分 异 时 间 和 演化 速率 ， 就 可 以 计算 在 给 定时 间 段 1 和 演化 速率 x 的 情况 下 ， 形 
态 特征 从 一 个 状态 变 为 男 一 个 状态 的 概率 ( 称 为 转移 概率 )。 这 个 概率 由 Mk 模型 (Lewis， 
2001) 给 出 。Mk 模 型 是 描述 特征 状态 变化 最 简单 的 模型 ， 它 假设 状态 之 间 转 换 的 速率 是 
相等 的 。 这 里 只 以 两 个 状态 的 特征 为 例 ， 用 Poo(r, 力 表示 状态 0 保持 不 变 的 概率 ，Poi(7, D 
表示 从 0 变 到 1 的 概率 ，Pio(7, 四 表 示 从 1 变 到 0 的 概率 ，Pii(7, 0 表示 状态 1 保持 不 变 的 概 
X, MPa, Ò= Py (r, À= 1/2 + 1/2 x exp(—2rt), P(r, À= P(r, À = 1/2 — 1/2 x exp(—2r1), 
从 公式 中 可 以 发 现 ， 时 间 t 和 和 速率 r 总 是 以 乘积 的 形式 出 现 ， 因 此 ， 在 没有 化 石 年 代 信息 
时 ， 两 者 是 不 可 识别 的 。 换 句 话 说 ， 单 单 依靠 形态 数据 来 建树 ， 树 的 校长 为 时 间 和 速率 
的 乘积 ， 即 距离 ， 以 每 个 特征 期 望 的 变化 次 数 为 单位 。 只 有 同时 利用 化 石 形态 和 年 代数 
据 才能 将 分 异 时 间 和 演化 速率 单独 估计 出 来 。 

以 图 2 为 例 ， 化 石 FI 和 F, 的 特征 状态 为 0, 现 生 类 群 S, 和 S, 的 特征 状态 为 1, 内 部 节点 的 
特征 状态 未 知 用 xo, xu x 表示 。FI 和 F, 的 时 间 分 别 为 100 Ma 和 50 Ma。 那 么 根据 Mk 模型 ， 
给 定时 间 树 7= {1, 世 和 速率 r 时 ， 特 征 状态 列 0011 的 概率 为 
P(0011|T, r) = Radel wall oe 站 Pet r)P. (ty r3) Poo (tits, r) Poi rs)P 0 (tits, r6) 

其 中 西格玛 符号 代表 对 特征 在 内 部 节点 所 有 可 能 状态 的 求 和 。 由 于 形态 特征 矩阵 
往往 只 包含 可 变 的 特征 ， 因 此 这 个 概率 还 要 除 以 所 有 可 变 状 态 的 概率 ， 即 P(0011 | 7, r)/ 
[1 — P(0000 | T, r)— PCL ID]。 带 有 这 一 校正 的 Mk 模型 称 为 Mkv 模 型 (Lewis, 2001). 

假设 形态 矩阵 中 的 特征 都 彼此 独立 ， 那 么 就 可 以 计算 每 一 列 特征 在 树 上 的 概率 ， 
再 把 这 些 概 率 乘 起 来 。 这 个 概率 被 称 为 似 然 函 数 ， 表 示 为 P(D |T, r), 其 中 D 代 表 形 态 特 
征 和 矩阵 数据 。 


5 贝 叶 斯 公式 


在 统计 推断 时 ， 参数 都 是 未 知 的 随机 变量 ,需要 根据 数据 来 估计 它们 的 分 布 ， 即 
计算 P(T, r, 9 | D), 该 分 布 称 为 后 验 分 布 ， 其 中 7 {r 世 为 时 间 树 ，r 为 演化 速率 ，0 代 表 
其 它 参数 (包括 4, ,wy 等 )。 根 据 分 层 贝 叶 斯 公式 ， 可 得 P(T, r, 8| D) = P(D| T, r) P(r) P(T| 4) 
P(0)/P(D)。 等 号 右 侧 分 子 中 ， 第 一 项 似 然 函 数 在 第 四 节 给 出 ， 第 二 项 演化 速率 的 先 验 分 
布 在 第 三 节 给 出 ， 第 三 项 和 第 四 项 为 时 间 树 及 其 参数 的 移 验 分 布 在 第 二 和 给 出 。 这 样 公 
式 分 子 中 各 项 都 可 以 计算 了 。 分母 P(D) 是 特征 数据 的 概率 ， 这 需要 计算 对 所 有 参数 的 多 
重 积 分 ,实际 上 基本 无 法 给 出 解析 表达 式 ， 只 能 通过 数值 算法 进行 近似 。 所 以 贝 叶 斯 计 
算 在 绝 大 多 数 情 况 下 会 使 用 马 氏 链 蒙 特 卡 罗 算 法 。 
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6 ” 马 氏 链 蒙 特 卡 罗 算 法 


马 氏 链 蒙 特 卡 罗 (MCMC) 算 法 通过 构造 马尔 科 夫 链 ， 使 其 平稳 分 布 为 要 估计 的 后 验 
分 布 。 这 里 为 了 简化 ， 只 以 一 维 参数 的 情形 为 示例 (图 3)。 实 际 分 析 中 ， 参 数 一 般 是 多 
维 的 (如 rz t, r, O), 不 过 算法 的 思想 是 类 似 的 。 

MCMC 采 用 的 Metropolis-Hastings 算 法 (Metropolis et al., 1953; Hastings, 1970) 可 分 为 
如 下 几 步 : 

a. 为 0 设 定 任意 初始 值 ; 

. 基于 0 当前 的 值 ， 建 议 一 个 新 的 值 0”, 例如 l ~ uniform(0 一 w/2, 0 + w/2); 
， 如 果 x(0”) > x(0), 就 接受 0”; 否则 ， 以 概率 a = aO) /x(0) 接受 0”; 

. WRC PES, UO = 0”; 否则， 保持 0 不 变 ; 

记录 6 的 值 ， 回 到 步骤 b。 

注意 到 ， 在 计算 概率 cx 的 时 候 ， 会 计算 后 验 分 布 的 比 ， 这 样 后 验 分 布 分 母 的 部 分 
就 约 择 了 ， 只 剩 分 子 部 分 的 比 。 也 就 是 说 ， 只 要 能 够 把 分 子 部 分 写 出 解析 表达 式 ， 
MCMC 算 法 就 可 以 使 用 来 估计 参数 的 后 验 分 布 了 。 

计算 结束 后 ， 就 收集 到 一 些 9 的 样本 。 由 于 参数 的 初始 状态 往往 比较 差 , MCMC 需 
要 经 过 很 多 代 才 收敛 到 后 验 概率 密度 比较 高 的 地 方 ， 因 此 在 估计 后 验 分 布 的 时 候 会 舍弃 
初始 的 一 些 样本 (burn-in), 只 用 MCMC 链 收敛 后 记录 的 样本 来 估计 后 验 分 布 。MrBayes 默 
认 会 弃 前 25% 的 样本 。 同 时 ，MCMC 链 还 要 迭代 足够 多 次 ， 以 保证 有 足够 多 的 有 效 样本 
来 估计 参数 的 后 验 分 布 。 一 般 需 要 有 效 样本 大 小 (ESS) 大 于 100。 

实际 运算 中 ， 最 好 独立 运行 至 少 两 次 MCMC, 以 确保 两 次 的 结果 是 一 致 的 。 有 时 
链 长 不 够 ， 或 者 不 同 的 运算 卡 在 不 同 的 后 验 分 布 区 域 ， 都 会 导致 估计 的 结果 不 一 致 。 
这 时 调整 MCMC 的 设置 或 者 改善 模型 都 可 能 帮助 MCMC 算 法 发 挥 更 好 的 效能 。 使 用 
Metropolis-coupled MCMC 也 是 有 效 跨越 多 峰 分 布 的 手段 (Lakner et al., 2008)。 该 算法 同 
时 运行 多 条 MCMC 链 ， 一 条 为 冷 链 (cold chain), 其 余 为 热 链 (hot chains), 热 链 和 冷 链 之 间 
可 以 相互 交换 。 当 然 MCMC 样 本 只 从 冷 链 中 采 


© S 


oo 


7O = 01D) 集 ， 热 链 只 是 用 来 帮助 跨越 多 峰 的 。MrBayes 
a(0) > a(O), accept 8» 默认 同时 独立 运行 两 次 ， 每 次 运算 使 用 四 条 
a(0) <O), l | 链 ， 其 中 一 条 为 冷 链 ， 其 余 三 条 为 热 链 。 
accept 0’ with [J | 

b.x(ONn(0) Pi: T 
= i y inital 7 讨论 
m : sh n 
0-w2 0 6+wi2 本 文 从 贝 叶 斯 统计 计算 的 角度 分 层 剖 析 
图 3 马 氏 链 蒙特 卡 罗 算 法 在 一 维 情 形 的 示例 。 了 支 端 定年 法 的 原理 和 计算 过 程 。 贝 叶 斯 后 验 


参数 9 的 后 验 分 布 z(0) 为 估计 的 目标 (曲线 ) ” 分 布 包 含 先 验 分 布 和 似 然 函 数 ， 其 中 先 验 分 布 
Fig. 3 Illustration of the MCMC algorithm for 的 两 个 重要 组 成 部 分 : 分 异 时 间 和 演化 速率 模 
one dimensional parameter Be R ` y Ż A ES E 
The posterior distribution x(0) of parameter 0 is 型 ， 在 定年 分 析 中 尤为 关键 ， 是 影响 定年 准确 
the target distribution to be estimated (curve) 性 的 主要 素 o 
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石化 生 灭 过 程 作为 描述 类 群 分 异 、 灭 绝 和 采样 的 随机 过 程 ， 具 有 较 大 的 灵活 性 。 
不 过 该 模型 还 有 竺 完善 之 处 。 现 生 类 群 的 采样 方式 可 以 是 随机 的 或 多 样 化 的 ， 这 两 种 情 
形 可 能 都 是 极端 ， 真 实 的 采样 模式 可 能 介 于 两 者 之 间 ， 或 者 有 的 支 系 是 随机 采样 ， 有 的 
文系 是 多 样 化 采样 ， 但 目前 还 没有 模型 能 够 支持 这 种 情况 。 在 更 好 的 模型 被 开发 出 来 之 
前 ， 可 能 只 能 调整 数据 来 尽量 符合 其 中 一 种 采样 策略 。 这 种 处 理 方式 在 现 生 类 群 比较 少 
或 根本 没有 时 一 般 问 题 不 大 。 男 外 ,分 异 速率 、 炎 绝 速率 和 化 石 采样 速率 可 以 按 分 段 的 
方式 随时 间 变 化 ， 不 过 在 同一 时 间 段 内 ， 所 有 树 梳 都 共享 同一 速率 。 对 于 不 同文 系 明 显 
受到 的 选择 压力 不 同 或 化 石 保存 的 完整 性 明显 不 同等 情况 ， 按 支 系 分 段 而 非 时 间 分 段 的 
模型 (Barido-Sottani et al., 2020) 可 能 更 合适 ， 虽 然 这 一 类 模型 本 身 也 有 其 他 限制 。 不 论 
怎样 ， 石 化 生 灭 过 程 只 是 作为 时 间 的 先 验 分 布 ， 当 数据 量 比较 大 时 (包括 化 石 在 树 上 的 
分 布 程度 和 数量 以 及 形态 特征 的 数量 和 完整 度 ), 数据 在 推断 中 会 起 主导 作用 而 先 验 的 影 
响 减少 。 但 是 实际 情况 往往 比较 复杂 ， 数 据 量 也 不 尽 如 人 意 ， 这 时 考察 不 同 先 验 的 影响 
就 尤为 重要 。 

演化 速率 的 先 验 ， 即 形态 钟 模型 ， 也 会 和 时 间 相 互 作用 ， 从 而 影响 对 最 终 分 异 时 
间 的 估计 。 这 种 影响 在 化 石 较 少 或 化 石 在 树 上 分 布 很 不 均 时 尤为 明显 。 这 主要 是 因为 
化 石 形态 数据 只 提供 了 距离 的 信息 (每 个 特征 期 望 的 变化 次 数 ), 其 为 时 间 和 速率 的 乘积 
( 见 第 四 节 )。 当 缺少 化 石 时 ， 就 没 办 法 准确 提供 时 间 的 信息 ， 那 么 对 于 同样 的 距离 ， 可 
以 是 很 长 时 间 速 率 很 慢 ， 也 可 以 是 很 短 时 间 但 速率 很 快 ， 具 体 是 怎样 就 只 能 取决 于 时 
间 和 演化 速率 的 先 验 了 。 对 于 某 些 支 系 时 间 估 计 得 明显 偏 大 或 偏 小 ， 但 又 没有 化 石 来 
校正 的 情况 ， 可 以 通过 添加 内 部 节点 的 校正 分 布 来 得 到 更 合理 的 时 间 佑 计 (O’Reilly and 
Donoghue, 2016)。 在 完全 没有 化 石 只 有 现 生 类 群 时 ， 节 点 定年 法 ( 另 一 类 型 定年 方法 ) 
(Yang and Rannala, 2006) 就 是 通过 使 用 内 部 节点 的 校正 分 布 来 完成 的 。 

描述 形态 特征 状态 变化 的 模型 也 有 很 大 的 改进 空间 ， 其 中 涉及 更 多 的 建 模 和 随机 
模拟 等 工作 ,不 是 本 研究 的 重点 ， 这 里 只 简单 讨论 一 下 Mk 模型 对 定年 可 能 的 影响 。Mk 
模型 假设 特征 各 个 状态 间 转 变 的 速率 都 是 相等 的 。 这 种 转变 有 无 序 和 有 序 之 分 (只 对 三 
个 及 以 上 状态 的 特征 )。 无 序 是 指 特征 可 以 直接 从 一 个 状态 变 为 任意 其 他 状态 (如 从 0 直 
接 变 为 3), 而 有 序 是 指 特征 只 能 在 临近 状态 间 直 接 变 化 (如 从 0 到 1, 从 1 到 2, 再 从 2 到 3)。 显 
然 ， 有 序 需 要 更 多 次 变化 (也 就 是 更 长 的 距离 ) 才 能 从 当前 状态 变 为 不 相 邻 的 状态 ， 因 此 
对 有 序 的 特征 假设 无 序 的 变化 会 低估 距离 。 更 复杂 的 情况 是 ， 各 个 状态 间 转 变 的 速率 未 
必 相 等 甚至 相差 很 多 ， 极 端 情况 像 Dollo 特 征 甚 至 是 不 可 逆 的 。 这 时 使 用 Mk 模型 也 会 造 
成 距离 估计 的 偏差 。 在 计算 时 还 假设 不 同 特征 之 间 都 是 独立 的 ， 如 果 有 些 特征 有 较 强 的 
相关 性 ， 则 会 导致 演化 距离 的 高 估 。 前 面 提 到 ， 上 距离 是 时 间 和 速率 的 乘积 ， 因 此 在 化 石 
(时 间 信 息 ) 很 丰富 的 理想 情况 下， 距离 的 偏差 会 主要 体现 到 演化 速率 上 而 对 分 异 时 间 影 
响 较 小 。 但 是 分 析 实 际 数据 时 会 更 复杂 一 些 ， 要 具体 问题 具体 分 析 。 相 关 的 工作 还 较 少 
(Klopfstein et al., 2019), 需要 更 多 后 续 人 研究 来 更 详细 地 考察 。 

最 后 提 到 贝 叶 斯 计算 使 用 的 MCMC 算 法 。 该 算法 的 策略 与 简约 法 和 似 然 法 有 明 
显 不 同 。 简 约法 寻找 的 是 简约 树 长 最 小 的 树 ， 似 然 法 寻求 的 是 似 然 值 最 大 时 参数 的 估 
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计 ( 即 最 大 似 然 树 )。 因 此 在 设计 树 的 搜索 方法 时 ， 只 要 尽 可 能 快速 地 找到 最 优 的 树 就 
可 以 了 。MCMC 算 法 是 为 了 估计 参数 的 后 验 分 布 ， 这 不 仅仅 是 一 个 点 ， 而 是 参数 的 空 
间 。 因 此 MCMC 算 法 的 效率 涉及 到 收敛 (convergence) 和 混合 (mixing) 两 部 分 。 收 敛 是 指 
MCMC 达 到 分 布 概率 密度 高 的 区 域 ， 混 合 是 指 MCMC 能 够 按 概 率 分 布 进行 取样 。 提 高 
收敛 速度 相对 容易 ， 可 以 通过 如 简约 树 长 向 导 的 方式 来 快速 找到 概率 大 的 树 (Zhang et 
al., 2020)。 提 高 混合 则 更 困难 ， 需 要 设计 更 好 的 建议 (proposal) 方 法 ， 这 是 贝 叶 斯 计算 的 
重点 也 是 难点 。 

总 之 ， 贝 叶 斯 支 端 定年 法 作为 整合 的 分 析 方 法 ,能够 结合 化 石 形态 和 年 代数 据 以 
及 现 生 类 和 群 的 形态 和 分 子 数 据 来 推断 类 和 群 的 系统 关系 ， 分 异 时 间 和 演化 速率 ， 同 时 考虑 
了 树 的 拓扑 结构 、 分 异 时 间 、 演 化 速率 以 及 化 石 年 代 的 不 确定 性 。 但 该 方法 仍 处 于 发 展 
初期 ， 模 型 和 算法 的 诸多 方面 还 吸 待 完善 ， 因 此 还 有 很 多 工作 需要 做 。 
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